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Abstract. We derive a physiologically structured multiscale model for biofilm development. The 
model has components on two spatial scales, which induce different time scales into the problem. The 
macroscopic behavior of the system is modeled using growth-induced flow in a domain with a moving 
boundary. Cell-level processes are incorporated into the model using a so-called physiologically 
structured variable to represent cell senescence, which in turn affects cell division and mortality. 
We present computational results for our models which shed light on modeling the combined role 
senescence and the biofilm state play in the defense strategy of bacteria. 
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1. Introduction. In this paper we derive a physiologically structured multiscale 
model for biofilm development. The model has components on two spatial scales, 
which induce different time scales into the problem. The macroscopic behavior of the 
system is modeled using growth-induced flow in a domain with a moving boundary, 
following m ^] . Cell-level processes are incorporated into the model using a so- 
called physiologically structured variable to represent cell senescence, which in turn 
affects cell division and mortality. We use "senescence" to mean "the organic process 
of growing older and showing the effects of increasing age" 1 . 

The multiscale nature of physiologically and spatially structured population mod- 
els, such as those in this paper and in differs from more typical multiscale 
systems where the smaller spatial scales have the faster time scales. In the structured 
multiscale systems, the dynamics of the relevant physiology of individuals within a 
population are homogenized to a distribution of a representative trait, such as age, 
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size or senescence. Although the underlying physiological system may have a very 
fast time scale (such as the protein network within a cell that controls the cell cycle), 
the distribution of the representative trait may evolve relatively slowly compared to 
the dynamics in space or in the reaction terms (such as the age distributions used 
to represent a tumor cell's position in the cell cycle in [2J, or a Proteus mirabilis 
multinuclear filament cell's length in pfl I12|). 

The derivation of the model in this paper follows that of Alpkvist and Klapper 
PP, with the addition of the explicit physiological structure in the bacteria popula- 
tions based on the notion of bacteria senescence demonstrated in [23] for cells with 
symmetric division and |SJ I2(J| for cells with asymmetric division. We also include 
explicit tracking of inert cell populations, which includes necrotic cells. 

To our knowledge, prior to this work, physiological structure has only been inte- 
grated into spatial models where motion is due to migration or taxis, represented by 
diffusion terms in the model equations |Hllfil ll2lllH| . Here, instead, motion is driven by 
growth- induced expansive stress, a much different mechanism that requires inclusion 
of a force balance equation. 

This paper is organized as follows. We first motivate the problem subject and 
derive the structured multiscale model of biofilm growth. Following this, we present 
a nondimcnsionalization and then a spatially homogeneous steady-state age distribu- 
tion, which in combination help to illustrate the differing age structures occurring 
in different places in the biofilm. Finally, we provide computational results for our 
models which shed light on the combined role senescence, and the self-organization 
into a biofilm state, play in the defensive capabilities of bacteria. 

2. Biofilms and Age Dependence. A biofilm is a collection of microorgan- 
isms, typically bacteria, enclosed within a self-secreted polymeric matrix. These films 
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are generally attached on one side to a solid boundary and, on the other, access sub- 
strates (e.g. oxygen) through a free surface, Figure ETT1 See [21J for a review. Biofilm 
properties change over long times (weeks) - it is possible to describe a maturation 
process in the development of a particular biofilm. That is, biofilms demonstrate 
aging effects. We are thus motivated to extend basic biofilm models to include age 
dependence. 

In fact, it seems that individual bacteria themselves suffer age dependence in the 
form of senescence. This property had been observed for some time in asymmetric 
dividers 8, 20j; more recently senescence has also been noted in the symmetric di- 
vider Escherichia coli |2H| . Under normal conditions, senescent cells make up a small 
percentage of the total population. However, aging (over medium time scales) may 
provide an effective defense against short time scale environmental disruptions but 
without affecting microbial community vitality during normal conditions. That is, 
we posit that multi-time scale behavior allows a powerful defense mechanism. In a 
recent paper |17j . it was argued that cell senescence offers a simple explanation for 
the phenomenon of bacterial persistence. Bacteria exhibit the phenomenon of "persis- 
tence" , the tendency for a small number of cells within a larger population to tolerate 
a wide range of antimicrobial challenges ^ n particular, the mechanism for 

this tolerance was suggested to be that senescent cells were less active and hence less 
susceptible than younger, more vigorous ones. Then, once the antimicrobial attack 
ceases, the persisting cells could produce new, vital cells which in turn would be ca- 
pable of regenerating the colony. Previously, others have argued that persisters were 
phenotypic variants [71 ITO! IT51 l2"Tl l2~2"] . 

Colony defense through persister cells is likely to be especially effective in biofilms 

where surviving persister cells, though perhaps small in number, have the opportunity 
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Fig. 3.1. Spatial domains for the senescence- structured biofilm model. 

to be protected by the biofilm matrix [T^]. As a result of this matrix they may 
have a particularly conducive environment for re-population once the antimicrobial 
challenge has ended. Thus, in order to demonstrate that senescent cells can distribute 
themselves throughout the biofilm and as a particular application of our age dependent 
biofilm model, we compute the spatial and temporal variation of age dependence. 
These senescent cells are produced on a medium time-scale (approximately 1 day): 
fast relative to the biofilm maturation time but slow compared to metabolic times. 
Thus persisters are generated quickly enough so that they can be an effective defense 
mechanism for mature biofilms but not so quickly that they interfere with competitive 
fitness. 

3. Derivation of the Model. We consider a spatial domain SI consisting of 
stratified subdomains B t for biomass and fl\B t for the bulk fluid. There are two 
moving interfaces in S7: T t separating B t from the rest of SI, and a bulk-substrate 
interface Tn b that is a fixed height iJf, above T t . The biofilm rests on a surface, 
denoted by a lower boundary, Tb- The spatial domains are illustrated in Figure mi 

We let bi(t,x, a), i — 1, . . . , iV&, denote the densities of the bacteria phenotypes 
in time t > 0, space x S B t , and senescence a > 0, and let denote their respective 
fluxes. The component of x representing height is denoted by z. Similarly, Cj(i, x), 
i = 1, . . . , N c , denote the substrate concentrations in time t > and space x € SI, and 



31 denotes their respective fluxes. In addition to active cell types, we allow for the 
presence of inert cells, including necrotic cells, that do not use or produce substrates, 
do not grow, and are not merely in a quiescent state. Lack of senescence allows us to 
list these separately from the bi because these cells do not have any age dependence. 
We let riiit, x) denote inert cells of type bi of all ages, and J™ their respective fluxes. 

Conservation of biomass yields equations 



dbi d(gi(a,ci,...,c N Jbj) 
dt da 



V • J? 



-fj,i(a,ci, . . . ,CAr c ,6i(t,x, •), . . . ,b Nh (t,x,-))bi(t,x, a) 

+ M&i(f,x,a), . . .,b Nb (t,x,a)), (3.1) 



for i = 1, . . . , Nb where /tj is the inactivation or "death" modulus with dependence 
on senescence, substrate concentrations and the densities of all bacteria phenotypes 
of all ages, and /, is the rate of net change to phenotypes i from all other phenotypes. 
The terms /; allows the possibility that bacteria have the capability of changing their 
phenotype in response to stimuli. A model that incorporated change due to mutation 
would do so in the age boundary condition. The senescence rate gi represents the 
physical wear-and-tear experienced by an aging individual in response to nutrient 
and/or oxygen uptake and exposure to waste. 

Due to the close relationship between senescence and chronological age, for this 
paper we make the simplifying assumption that = 1, for i — 1, . . . , Nb. Further 
below, we will specify senescence as a function of age for use in both the inactivation 
modulus, jl, and the fecundity, 0, defined below. Equations (|3.1f) then become the 
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age- and space-structured equations 

§ + ^ + v.j^ 

at oa 

= -fj,i(a,ci, . . .,c Nc ,bi(t,x,-), . . .,b Nb (t,x, ■))&*(*, x, a) 

+ Mbi(t,x,a),...,b Nb (t,x,a)). (3.2) 

Observations in |23| showed that even in symmetric cell division, one of the two 
new cells contains older material and overall less vitality than the other (referred to 
as "old pole" and "new pole" cells, resp.). This results in a physiologically struc- 
tured mathematical representation of the bacterial cell cycle that is closer to that 
for birth-death processes in animals than what has been commonly used to repre- 
sent cell division [HHJS]. These models were built on the assumption that a mother 
cell divided into two daughter cells of equal and high vitality, thereby assigning each 
daughter cell a senescence of zero and removing the mother cell from the population. 
In our old-pole/new-pole formulation rather, the mother cell remains in the popula- 
tion and continues to undergo senescence from the point it had at cell division while 
giving rise to a single daughter cell with senescence zero. This notion of senescence 
allows two implications, based on the account in |23| . that are relevant to the model in 
this paper. First, old-pole cells grow slower than new-pole cells produced in the same 
division. Second, old-pole cells become inert at a higher rate than new-pole cells. 

Although a more elaborate model would include explicit size structure similar 
to what was done in models in [S] 125) . we can take advantage of continuous senes- 
cence and time to incorporate the first old-pole/new-pole issue into a fecundity term, 
$i(a, ci, Cjv , b±(t, x, •),..., 6ATi,(i,x, •)), with dependence on age, substrate con- 
centrations and the densities of all bacteria phenotypes of all ages. Differences in 

individual sizes influence volume fractions, represented by giving mother cells with 
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larger offspring a corresponding higher fecundity. The fecundities, 0, account not 
just for differences in daughter size due to the mother's size, but also heterogeneities 
in the mean growth rates across phenotypes i. A third more minor property of cell 
senescence mentioned in |23| , that new-pole cells are marginally more likely to divide 
sooner than old-pole cells, can also be included in (3. (Similarly, we can incorporate 
the second property of higher incidence of becoming inert into The resulting 

senescence boundary condition, or "birth" condition 2 , is 

/•OO 

&i(i,x, 0) = / j3i(a, . . .)bi(t,x,a) da, for i = 1, . . . , JVj. (3.3) 
Jo 

We represent retention of inert cells using Nb inert-cell classes governed by the 

conservation equations 

dn (t xl f°° 

iV ' ; fV-Jf= fri(a,...)bi(t,xL,a) da, for i = 1, . . . ,N b . (3.4) 



Conservation of substrate mass yields 

^+V-J? = r i , for j = l,...,iV c , (3.5) 

where rj denotes gain or loss of the j-th substrate concentration through interactions 
with the biomass such as consumption or excretion. Assuming Fick's Law gives J'j = 
—DjVcj for constants Dj. The substrate masses are also subject to advection, but 
the velocity is sufficiently slow that we can neglect the advective contribution to the 
flux. Likewise, substrate material diffuses several orders of magnitude faster than 
the rates at which bacteria grow or advect, allowing us to make a quasi-steady-state 
assumption so that 

-DjV 2 c 3 = rj , for j = 1, . . . , N c . (3.6) 



2 A model where transition between some of the different classes occurs due to mutation would 
have an age boundary condition of the form bi (t , x, 0) = Yl k ^ik Jq° Pk^k (*> x i a ) where M ik is a 
matrix with mutation rates in the entries not on the main diagonal, and one minus the sum of those 
rates on the main diagonal. The specifics of the entries of Mik would account for what mutations 
underlie the different phenotypes. A version with linear progression through phenotype classes was 
used in UJ. 



We let x, a) and Pi(t, x, a) denote the volume fraction per age and density per 
age relative to volume fraction, resp., of phenotypes i, so that bi = piOi. We assume 
incompressibility of biomass with pi(t,x) = p* for positive constants p*. We also 
assume inert cells have the same incompressibility properties, and the same densities 
relative to volume fractions, p*, as active cells. We let rji(t,x) denote the volume 
fraction of inert phenotype i cells, which is related to the density of inert phenotype i 
cells by rii = p\r\i. We assume such cells all behave the same regardless of phenotype, 
and track only the total volume fraction of inert cells of all phenotypes, denoted by 
Af(t, x) . Equations (|3.4|) , rewritten as 

d Vi (t x) + J_ v r°° fori=l,...,JV 6 , (3.7) 

ot Pt Jo 

become, after summing over i, the governing equation for Af, 

2 _V-J? = J M(*,x), (3.8) 



rW((,x) " 



<9f 



where 



^6 

AA(t,x) =^^(t,x), (3.9) 

■^t /.CXD 

M(t,x)=y2 p l {a,...)^ l (t,x,a) da. (3.10) 



i=0 

We require the biomass volume fractions to total to one so that 

■Ni> poo 

7V(t,x) + V/ i?i(t,x,o) da= 1. (3.11) 

i=i ^ 

Assuming that transport of biomass, including inert cells, is governed by an ad- 
vective process, with a volumetric flow u(t, x) for all classes and ages, gives the fluxes 
3\ = p*i?jU for t=l..., Nb- Following we assume that the volumetric flow is 

stress driven according to 

u = -AVp, (3.12) 



where p(t, x) is the pressure and A > the Darcy constant. As in P|^, p = 
in fl\B t . Pressure is determined in order to enforce incompressibility in response to 
growth (see below) and hence (|3.12|) can be viewed as a balance of growth-induced 
stress against friction. Other choices of force balance are possible. 

Substituting bi — p*di and J\ = p*$iU into equations 1|3.2|) gives, for i = 
l,...,N b , 

— + + V ■ (ui?j) = -m{a,cx, . . . ,c Wc ,i?i(t,x, •), . . . ,<& Nb (t,x, ■))i3 l (t,x,a) 

+ / i ( 1 ? l (i,x,a),...,^ 6 (i,x,a)), (3.13) 

where 

Hi(a,ci,..., c Nc , ^(t, x, •),... , i} Nb (t, x, •)) = 

Ai(a,ci, . . .,c Nc ,bi(t,x, •), ■ ■ ■ ,b Nb (t,x, •)) (3.14) 

and 

/i(#i(i,x,a), . . . , d N (t, x, a)) = ^-/i(&i(t,x, a), . . . , x, a)). (3.15) 
The birth conditions (|3.3|l become 

/>OG 

i?i(t,x,0)= / &(«,... )0i(t,x, a) da, for * = 1, . . . ,JV 6 . (3.16) 
Jo 

where 

A(<I,Ci, . . .,CAT o! 1?l(t,X,-), ■ ■ ■ ,$Ni(t,X, •)) = 

/3i(a,ci, . . . ,c Ne ,bx(t,x,-), . . . ,b Nb (t,x, •)). (3.17) 

Substituting J™ = p*rjiU into equation (|3.8(l and using equation (|3.9(l gives 

^— =A4-V-uAf, (3.18) 
at 
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Integrating equation (|3.13|) over age and summing over i gives 

d (sr r a a wV/r / ,00 a? i (t,x,a) 



da 



=v-(uJV)-A< as denned below 

■ V- ( y $i(i,x,a) da J - y ^(a, . . . 7. x. a i </a 



=V-(u(x-JV0) =A1 

,i9 Nb (t : x,a)) da I . (3.19) 



/ ^6 /-OO 

+ V/ /i(i?i(t,x,a),. 



=;f as denned below 

Using equations (|3.11(l and H3.18[l . we find that the first term in the first line of 
equation 1)3. 19|) is dt(l — TV) = V • (uN) — M.. For the second term in the first line, we 
assume that for i = 1, . . . ,Nb, are sufficiently smooth, and that the corresponding 
fii are bounded away from zero for a large, so that each di will eventually decay 
exponentially to zero as a — > oo (see section 7 in 0). We then obtain for the negative 
of the second term, using the age boundary conditions defined by equation (|3.16(l , 

B(i,x)=^^(t,x,0) 

i=Q 

= J2 Pi{a,...)di(t,x,a) da. (3.20) 
i=o Jo 

For the first term in the second line of equation l|3.19[l . we again use equation l|3.11|l 
to obtain V • (u(l — AT)). Recall that the second to last term of equation (|3.19() is just 
M(t,x) and set the last term to 

N b poo 

^(*,x)=V/ fida. (3.21) 
i=a Jo 

We note that T is not generally identically zero; a similar sum over all phenotypes 

of the integrals over all ages of the net changes between phenotypes, /, is conserved 

to be zero. However, because the densities relative to volume fractions, p*, are not 

identical, we generally only have T = when p* = p* for some constant p* and for 
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alH = 1, . . . , Nj,. We rewrite equation l|3.19fl more compactly as 

V ■u = B(t,x)+T(t,x), (3.22) 

the incomprcssibility relation for our system. 

Substituting u = — AVp results in an equation for the pressure, namely 

-\V 2 p = B(t,x) + T(t,x), mB t . (3.23) 

Distributing the divergence operator, and again using u = — AVp along with equation 
l(S23l, gives us V ■ (ui?j) = -AVp ■ + + J 7 ), so that equation lj3~T3|) can be 
rewritten, for i = 1, . . . , Nf,, 



d-d di) 

_^ + —L _ AV p ■ Wi = + /„• - + n (3.24) 
ot aa 

Similarly, we rewrite equation l|3.18|l as 

= M + XVp-VAf -N{B + T). (3.25) 

We see from equation (|3.23|) that p is proportional to A -1 , so that AVp is independent 
of A. Consequently, di and N are independent of A, allowing us to set A = 1. 

We impose periodic and other boundary conditions similar to what was done in 
P to obtain the complete model, for i = 1, . . . , Nb and j — 1, . . . , N c , 
fid- 

-^ 1 + ^-Vp-Vi? i = -^+/ i -i? i (B + ^), x£B t ,i>0,a>0, (3.26a) 
ot Oa 



poo 

/ Pi(a, . . .)-di(t,x,a) da, xeB t ,(>0 
Jo 



$i(t,x,0) = I /3i(a,...)di(t,x,a) da, x€B t ,t>0, (3.26b) 

d-d- 

—^ = 0, x e T B ,t > 0,a > 0, (3.26c) 

Oz 

i? i (0,x,o) = ^(x ) o), xeB t ,a>0, (3.26d) 

^-Vp-VAA = X-AA(S + Jc), xe5,,i>0, (3.26e) 

OT 

^-=0, xer B ,t>0, (3.26f) 
Oz 

M(0,x) = N°(x), xeB t , (3.26g) 
11 



-V 2 p = B + T, x6B,,t>0, (3.26h) 
P = 0, xer f! (>0 ; (3.26i) 
^=0, xer B ,(>0, (3.26j) 



-DjV 2 Cj =rj, x e n,t > 0, (3.26k) 

JV=0, xefl\B ( , (3.261) 

cj = c*, xer ffi ,(>0, (3.26m) 
<9c- 

-7^=0, xer B) (>0. (3.26n) 

oz 

The normal velocity of the interface is given by 

-Vp-n=-|£ (3.26o) 

where n is the unit outward normal of Tb ■ 

We make particular choices of functions /3 and /i as follows. To reflect the dimin- 
ished new-cell production by senescent cells discussed in , we define senescence as 
a function of age, <r(a), such that <r(0) = and a (a) — > 1 as a — > oo, and incorporate 
cr(a) into /3(a, c) and /z(a, c), 

/3(a,c)=A)(c)(l-a(a)), (3.27a) 
/i(a, c) = jHo(c)(r(a). (3.27b) 

We neglect the c dependence of jtxo, and choose 

o-(o) = -^— , a>0, (3.28a) 
a* + a 

A>(c) - ir~r~ ' ( 3 - 28b ) 

where a* is the senescence age scale, and, following [Q, is the maximum growth 
rate (with units of inverse age) and k is the Monod saturation constant (with units 
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of concentration). Oxygen uptake has the form 



r(c,i?(t,x, •)) = 7^- / (l-a(a))&(t,x,a) da, (3.29) 



k + c Jo 

where £ is the maximum uptake rate. 

4. Non-dimensionalization. We simplify in the following to Nf, = N c = 1, i.e., 
restrict to one active phenotype and one substrate, and drop indexing subscripts. Note 
now that N — r](t, x) and M — J °° [i(a, ■ ■ .)${t, x, a) da. Also note that / = T = 0. 
We will continue to assume that /3 = /3(a, c) and = /i(a, c). 

Let /3 be a typical value of /3(a, c) and let p, be a typical value of /i(a, c). Choosing 
a characteristic time scale T = 1//3, age scale A = 1//2, and temporarily reintroducing 
the friction coefficient A, we nondimensionalize according to t — t/T, a = a/ A, x — 
x/L, and = f3T, p, = pA, c = c/c*, r = r/r(c*), p = p{XT/L 2 ), = fj = ??. 
Here L is a (problem-dependent) characteristic system length scale. 

Substituting into the system l|3.26(l and dropping tildes, we obtain 

dti dd f°° 

— +A— -Vp-Vtf = -A/z0-tf / (3(a,c)ti{a)da, x e B t ,t > 0,a > 0, (4.1a) 
at oa Jq 



/>oc 

i?(i,x,Q)=/ /3(a,c)*d(t,x,a) da, xeB t ,i>0, 
Jo 



(4.1b) 



9?? 

— = 0, x e r B ,t > o,a > o, (4.ic) 

i?(0,x,a) = i?o(x,a), xeBi,a>0, (4.1d) 



a /-oo poo 

- Vp- V77 = A / u(a,c)'&(a)da- r> f 0(a,c)-3(a) da, xeB t ,i>0, (4.1e) 



5r = o, xer B ,t>o, 

j?(0,x) = ??o(x), x€B(, 
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(4.1f) 
(4-lg) 



V 2 p=-/ P(a,c)-d(a)da, xeB t ,t>0, (4.1h) 
Jo 

P = o, xer ( ,i>o, (4.ii) 

^ = o, xer s ,t>o, (4.ij) 
az 



V 2 c=-Gr, xe!],i>0, (4.1k) 

r = 0, x€f!\B t , (4.11) 

c=l, xer flt ,i>0, (4.1m) 

-^=0, xer B ,i>0. (4.1n) 



Here G = L 2 r(c*)/(c* D) and thus 1/vG, the active layer depth, is a non-dimensional 
measure of the depth (scaled by system size) to which substrate can penetrate into 
the biofilm before it is consumed. Likewise 

A = i (4.2) 

is a non-dimensional ratio of characteristic deactivity time to characteristic reproduc- 
tion time. 

The non-dimensional forms of equations (|3.28a|l and l|3.28bjl are 

a(a) = — , a > , (4.3a) 
b + a 

A(c) = (4.3b) 
K + c 

where S — a* /A is a comparison of senescence age with system age scale, K — k/c* is 
a measure of saturation level (large K means substrate-limited behavior and small K 
indicates growth- limited behavior), and P — tp/ f3 is a measure of maximum to typical 
yield. 

The magnitude of A may depend on location within the biofilm. We identify two 

regimes. First, near the top of the biofilm, in particular within the active layer, c is 
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0(1) and we can then generally expect for a viable biofilm that (3 be large compared to 
Jx, i.e., A small. In this case advective terms dominate in equations (14. laf) and H4.1c(l 
over the death terms. If advection is unimportant, i.e., Vp • V$ is small, then age 
scale is determined by the second term of (|4. 1 al) which then requires d/da ~ A -1 . 
Such scaling is in fact observed within the biofilm active layer, see Section 

Second, beneath the active layer, equation l)4.1k(l indicates exponential decay (in 
space) of c. Hence, below a sharp transition region from the active layer, we can 
expect (3 to be small compared to ft, i.e., A large. In this case at first glance equation 
l|4.1a|l indicates i) has an /i-governed decaying age structure. There is a subtlety here 
however. Large A in (|4. suggests an approximately exponential age profile of the 
form f?(a) = i?(0) exp(— /ia). But such a form is inconsistent with l|4.1bj) . which does 
not allow dependence on /i, unless $(0) = 0. In other words, for large A the birth 
term is insufficient to introduce enough new cells to overcome death and so an active 
population is not viable. Having said this, however, we will observe a /i-determined 
exponential age structure develop in the deeper parts of the biofilm, see Section El 
The reason for this is that in the lower layer of the biofilm, where the population is 
barely viable, the birth rate has decreased to the point that it is only just balancing 
death. Hence condition 14.1b|) . which requires that an exponential age structure be 
controlled by /3, also implies that exponential age structure be determined by [i. 



5. Spatially Homogeneous Steady-State Age Distributions. We assume 

spatial homogeneity and temporal stationarity, i.e., •& — f?(a), r\ — 77(a), on —00 < 

z < 00. (We note that pressure gradients within the biofilm and hence advection are 
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generally weak within inactive regions.) Then rj satisfy 



dfi 

— = -^ + $ )§ 1 (5.1) 
da 



oc 



d a = / (3d da, (5.2) 
Jo 

/>OC 

rj = 1 — / rfi, 



where t?o = i9(0). This system is unphysical in that it requires unbounded velocities 
(the pressure gradient takes the formp z = C1Z + C2) to enforce incompressibility, but 
it is nevertheless useful for illustrative purposes. 
The solution to equation 1|5.1|) is 

0(a) = ti Q e-5> da 'e- &oa . (5.3) 

Thus equation l|5.2(l implies the condition 

/>oo 

190=/ /3tV~ io<VdQ V' ,oa da. (5.4) 
Jo 

In order to have a nontrivial solution, we require $0 to satisfy 

/•OO 

1= / (3e~foH da ' e -»oa da ^ ( 55 ) 



II 



if possible. If this is not possible, then $0 = i)(a) = is the only solution. 

The choice of /i and (3 independent of a allows a particular transparence. In this 
case condition (|5.5|l becomes 



1= / pe-(» + ° o)a da, (5.6) 
Jo 

which has a solution with tDq > if 

/ (3e- (ia da>l, (5.7) 
Jo 

i.e., if new cells can be produced sufficiently fast to replace aging (and dying) ones. 

Note that this condition cannot be satisfied for // sufficiently large or (3 sufficiently 
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small (in which case $0 = necessarily). Now if we write $0 = i9 Q — fx, then $ solves 



1 = 



/ (3e-® aa da, 



(5.8) 



Jo 



that is, $0 = 0- Equation (|5.3|l becomes 



d(a) = {0-^ a . 



(5.9) 



Fecundity fixes the profile of the age structure though age distribution amplitude 
depends on both (3 and /i. Large results in a steep age profile with amplitude 
almost independent of /1. Small results in a flat age profile. However we note that 
the viability boundary (in parameter space) occurs at = \i. For marginally viable 
populations — \i + e and hence the population age profile is exponential with decay 
rate approximately /i. Note that n gives the slowest possible rate of decay. With 
regards to a biofilm model, if we assume that decreases with decreasing c, then age 
structure should flatten deeper down into the biofilm. In fact, we expect an abrupt 
transition from steep to flat profile as we pass through the active layer. 
Returning to our specified forms of 0(a, c) and fi(a, c) we have 



Going back for a moment to dimensional variables, we use days as units of time, take 
/io = 0.25, and vary 0o(c) to induce changes in a characteristic reproduction rate, 
= la ^( a ' c ) W e set a* = 0.5, which accounts for a loss of 1% vitality per 
division (occurs on average every 0.02 days for E. coli). Figure 1ST! shows the response 
of the steady-state solutions to changes in 0. We note that the steady-state active-cell 
distribution is zero when is less than roughly 0.33. 




(5.10b) 



(5.10a) 
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1 2 3 4 U 1 2 3 

age (days) Characteristic p 

(a) Response of ln(0(a)) to changes in j3. (b) Response of r\ to 

changes in /3. 

Fig. 5.1. Response of steady- state solutions to changes in the characteristic reproduction rate, 
j3 = Jq j3{a, c) da. 



6. Computational Results. In this section we present computational results 
for one spatial dimension (height of the biofilm) , and explicit age structure represent- 
ing cell senescence, for the dimensional system (|3.26|l . The height of the biofilm, T t , is 
regulated using an erosion term at the biofilm/substrate interface (a standard devise 
in biofilm models, see e.g. |15|) given by 



dT t _ dp 
dt dz 



-aTi, (6.1) 

z=r t 



where a is the erosion coefficient. 

For the computations presented in this section, we consider again the case of 
TVf, = N c = 1. We take as the initial condition a biofilm with a height of T t (0) = 50/im 
and with an age distribution that is initially the same for all heights, i?(0, z, a) = 
0.35 * max(l — ^, 0) for < z < 50/^m. This piecewise linear function, when converted 
from age structure to senescence structure (recall cr(a) — a «°^ a ) 3 , closely approximates 



3 Computations with different a{a), namely a{a) = 1 — exp(— a/a*), c(a) = max(j^j-, 0.999) 
and cr(a) = max(^, 0.999), and with the same parameters as in this section, except for if), give 
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Time=0 days Time=5 days Time=10days 




Fig. 6.1. Biofilm dynamics from initial colonization to steady state. The height of the colored 
area represents the height of the biofilm, including both active and inert bacteria. Color represents 
cell state: black represents inert cells, and a spectrum from off-white to yellow to orange to red 
represents senescence of a cell of a given age, a(a). The horizontal width of a color constitutes the 
volume fraction of cells of the corresponding senescence. 



the senescence structure at the top of the biofilm when it is near steady state and 
thus represents a situation where a new area is being colonized by material from the 
top of a mature biofilm when it is near steady state. The motivation is that a young 
biofilm may be formed by colonization of cells detached from the upper region of an 



upstream, mature biofilm. 



We use a time unit equal to one day. As a result, we take the senescence time 
scale to be a* = 0.5, as was done in Section [SJ We use the division time for E. coli, 



which is roughly every 30 minutes. 



We set the erosion parameter to be a = 0.03, the distance between T t and Tjj b to 
be Hb — 37.5/im, and the parameters for the various functional forms to be fio — 0.25, 
k — 0.05, ip = 2, £ = 3, and c* = 1. These parameter values are of the same order of 

qualitatively similar results. We need to change if) since the different areas under the curves of <r(a) 
give substantially different total new cell production. 
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age (days) 

Fig. 6.2. Normalized age distributions, ignoring the inert cell populations, one third of the 
way from the bottom, two thirds of the way from the bottom, and at the top of the biofilm at time 
t = 25 days. The plot o/#(a) = 0.01 exp(— a/4)(l + 2a) 1 / 4 is the large A limit of equation l5.1Uat . 
The coefficient of 0.01 governing the magnitude of the curve is chosen for ease of comparison. The 
plot o/i9(a) = 0.7676 cxp(— 0.7471 a)(l + 2a) 1 / 8 is the re-normalized steady-state, given an oxygen 
concentration of c = 0.5581, in the absence of advection. Differences between this curve and the 
computed solution at the top of the biofilm illustrate the role of advection, including the upward flow 
of a relatively greater proportion of inert and senescent cells 

magnitude of those used in pQ, with modifications due to the inclusion of age structure 
in the model equations. 

Results of the computations with the above parameters are displayed in Figure 
16.11 The height of the colored area indicates the height of the biofilm, including both 
active and inert bacteria. Color represents cell state: black designates inert cells, and 
a spectrum from off-white to yellow to orange to red designates senescence of a cell 
of a given age, a (a). The horizontal width occupied by a color indicates the volume 
fraction of cells of the corresponding senescence. 

The biofilm tends to a steady state, as discussed in Section [5J consisting of an 
active layer at the top and passive layer appearing as a stalk. It is already understood 
that this physical structure provides a form of protection for the bacteria population 
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as a whole .9 ■ The question remains: how does senescence and the corresponding 
resistance to antimicrobial challenge fit into the overall defensive strategy of bacteria? 

To obtain an answer, we first consider the normalized age distributions, ignoring 
the inert cell populations, one third of the way from the bottom, two thirds of the 
way from the bottom, and at the top of the biofilm at time t = 25 days. These 
distributions are shown in Figure E2I 

As we descend the biofilm down through the active layer and into the passive layer, 
we expect A to increase so that equation <|4.1a|) approaches, in steady state, equation 
I5.10al The plot of d(a) = 0.01 exp(-a/4)(l + 2a) 1/4 highlights the convergence 
toward the shape of the curve of the large A limit of equation l|5.10a|l . The coefficient 
of 0.01 governing the magnitude of the curve is chosen for ease of comparison. At 
the top of the biofilm at steady-state, the oxygen concentration is approximately 
c = 0.5581. Using this value, our specific functional forms defined in equations (|3.27l) - 
H3.28fl . and equations l|5.3|) and (|5.5|l . we obtain a value of 9q — 0.6221. We re- 
normalize the age distribution to total one so that equation (|5.3|) has the specific form 
■#(a) ~ 0.7676 exp(— 0.7471 a) (1 + 2a) 1 / 8 . This represents the situation when there 
is no advection. Differences between this function and the graph of the computed 
steady-state age distribution at the top of the biofilm reflect the role of advection, 
including the upward flow of material with a relatively higher proportion of inert and 
senescent cells. 

We find the expected result, as discussed in Sections 01 and [S] that the age dis- 
tributions broaden as we go from the active to the passive layers within the biofilm. 
In the inactive region, the profile matches that of an approximately growth-death 
balanced population. This is to be expected in an erosion maintained steady state 
- some growth must occur all the way to the bottom of the biofilm. A non-viable 
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zone does not form in the presence of erosion because such a zone does not result 
in any growth induced pressure and hence does not increase expansive velocity. We 
remark that it is possible that /i-dominated age structure in the inactive region of the 
biofilm is thus a byproduct of erosion. A more general (and realistic) biofilm model 
would allow for the possibility of mechanical detachment; this would require a much 
more elaborate set-up than the one used here (i.e., mechanical stress coupling in three 
dimensions). It is still plausible that a non- viable zone would lead to detachment and 
so we posit that age structure would not change. 

Finally, we extend the computation to include the effects of antimicrobial chal- 
lenge. We assume the antimicrobial agent has a source at the bulk-substrate interface 
Th,, ) that it diffuses on a fast time scale compared to growth, and, for simplicity, 
that it is not degraded by the biofilm. Consequently, the antimicrobial saturates the 
biofilm essentially instantaneously and thus we can model the effects of antimicro- 
bial challenge by modifying the death modulus /x rather than adding an additional 
chemical species equation to our system: 



This form assumes that older cells are more resistant to antimicrobial challenge than 
younger cells, and that the antimicrobial agent affects metabolically active cells more 
than less active cells, represented by the oxygen dependence of fi\. In particular, we 
take 



Results are shown in Figure RHfl for the case when an antimicrobial agent is applied 
from time t — 35 to time t = 35.2. The senescence structure of the population in 
the stalk allows it to maintain itself even after the active layer is largely decimated. 



fj,(t, a, c) = (J, a(a) + /ii(t, c)(l - a (a)). 



(6.2) 




35 < t < 35.2, 
otherwise. 



(6.3) 
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Time=35 days Time=35.2 days Time=35 4 days 




0.5 10 0.5 10 0.5 1 

cumulaiive volume fraction cumulative volume traction cumulative volume traction 



Fig. 6.3. Response to an antimicrobial agent applied from time t = 35 to time t = 35.2. The 
height of the colored area represents the height of the biofilm, including both active and inert bacteria. 
Color represents cell state: black represents inert cells, and a spectrum from off-white to yellow to 
orange to red represents senescence of a cell of a given age, a(a). The horizontal width of a color 
constitutes the volume fraction of cells of the corresponding senescence. 



Moreover, upon removal of the antimicrobial agent, the population of older cells in 
the active layer is quickly replaced by younger cells, which then return the biofilm 
to its steady state over a longer maturation time. Note that the height continues to 
drop after removal of the antimicrobial agent prior to regrowth. 



6.1. Numerical Methods. We employ a moving-grid Galerkin method in age, 
using piecewise constants as the approximation space [2]. The use of higher-order 
approximation spaces in age was discussed in 0]. The moving-grid Galerkin method 
decouples the age and time discretizations, while allowing age and time to advance 
together along characteristic lines. Consequently, we are able to solve the model equa- 
tions in age without numerical dispersal or oscillations. Because the transport in age 
is computed by the movement of the grid, rather than by a difference approximation 



of the age derivative or through jump terms in a standard discontinuous Galerkin 
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method, the only meaningful source of error is approximation error, which underlies 
the superconvergence results in [2 El- For the case of piecewise constant functions, 
we obtain a second-order correct method in age. 

We integrate time using a step-doubling method 5 . Step-doubling consists of 
taking one step of backward Eulcr over a time step, and then taking two half steps of 
backward Euler over the same time interval. This results in two things. First, we can 
compare the two late-time solutions for the error control needed for the adaptivity 
in time. Second, we can extrapolate the two solutions to get a likely second-order 
accurate solution in time. 

For the spatial variable, we discretize, over a uniform partition, the domain [0, Ft] 
and compute the changes in biofilm height by solving equation (|6.1|) . We impose a 
boundary condition on c at T t by using a ghost node positioned at Tu b - Fluxes are 
computed using upwind differencing |18| . Although this method is only first-order 
correct, it has been sufficient for the computations presented in this paper, given the 
lack of sharp fronts in the interior of the biofilm, [0, T t ]. More advanced methods will 
be needed for computations with more spatial dimensions. 

The discretizations in the computational results presented above used a uniform 
partition of the spatial interval [0, T t ] with 301 nodes, and a uniform age discretization 
of the truncated age domain, [0,16], with Aa =1/8 and piecewise constant basis 
functions. A uniform age discretization in the context of the moving-grid Galerkin 
method means that all but the first and last age intervals are constant in length, and 
that a new age interval is introduced at the birth boundary when the old birth interval 
reaches Aa in length. The tolerance parameter for the adaptive time-stepping in the 
step-doubling algorithm was 5 x 10~ 3 . 
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7. Conclusions. In this paper we presented a multiscale model and simulation 
of biofilm development that is interesting for three major reasons. One is the nonstan- 
dard multiscale nature of the problem: cell division and aging is a result of complex, 
and fast, micro- and nano-scale processes, at least when compared to the advective 
scale of the biofilm growth. However, by representing the cell division and aging 
process using notions of senescence and age, we have a mechanism for the cellular 
scale that, in keeping with what has been observed in other age- and space-structured 
multiscale systems El El El, is in general slower than the advective process. But 
even here we see novelty; unlike El El Ej j the relative ranking of the time scales 
of the aging and advective processes inverts as we move from an active layer at the 
top of the biofilm to a passive layer below. Further, both of these times scales are 
fast with respect to the biofilm maturation time. 

This inversion of the time scales underlies another major point of interest: the 
implication that the active layer does not merely provide a physical shield for a reser- 
voir of cells in the passive layer, but also induces the passive layer to consist of an 
increased proportion of senescent persister cells. 

A third point of interest, and one which may have relevance to other biological 
systems that exist in a polymer matrix, e.g. tumor-matrix interactions, is the novel 
inclusion of age structure in a spatial model where movement is due to growth-driven 
expansive stress rather than diffusion or diffusion-like terms that represent mecha- 
nisms such as chemotaxis or haptotaxis (movement of cells up a matrix gradient). 

The model in this paper has a number of entailments for future work. One is 
experimental verification of the hypothesis that passive layers in biofilm contain a 
disproportionate number of persister cells. Another is a generalization of the model 
to higher spatial dimensions, and a study to see in what manner the physical stalk of 
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the mushroom-like shapes biofilm often form affects the persister "stalk" visualized in 
the senescence structure in this paper. Finally, it is likely that many of the modeling 
and simulation ideas developed in this paper have relevance to other systems. For 
example, inclusion of growth-driven expansive stress into an age- and space-structured 
tumor model like that in alongside other major mechanisms of motion such as 
diffusion and haptotaxis, would result in models with more fidelity to the physical 
mechanisms of tumor invasion. 
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